#################################################
# Section 5: Descriptive statistics
################################################
out_path<-paste0(wkdir,"/output_uncorr/2_descriptives/")
if(!dir.exists(out_path)){
	dir.create(out_path,recursive=T)
}else{
	warning("Overwriting previous run, if any")
}

############
#Descriptive stats for variables:
#
#Middle class preferences
#Affluent preferences
#Business interest group preferences
#Advocacy interest group preferences
#Foreign, economic, and social policies
#Ideological codes of policies
#Distribution of policy passage
#Republican support for policies
#Democrats support for policies
############
setwd(out_path)
vars<-aic_clean[,c("pred50_sw","pred90_sw","b.scaled.intgrp","m.scaled.intgrp",
					"fp.policies","econ.policies","ne.policies","rescaled_ideol_code","bioutcome",
					"rep","dem")]
					
desc.table<-psych::describe(vars,check=T)
rownames(desc.table)<-c("Middle Class opinion","Affluent opinion","Business int grp opinion","Advocacy int grp opinion","Foreign policies ({0,1})","Economic policies ({0,1})","Social policies ({0,1})","Ideological code","Policy passage ({0,1})","Republican support","Democratic support")
desc.table$trimmed<-NULL
desc.table$mad<-NULL
desc.table$se<-NULL
desc.table$skew<-NULL

sink("Descriptives.tex")
print(xtable(desc.table))
sink()

###########
#Correlation
###########
.z<-aic_clean[,c("pred50_sw","pred90_sw","b.scaled.intgrp","m.scaled.intgrp",
					"rescaled_ideol_code","bioutcome","rep","dem")]
.z<-.z[complete.cases(.z),]
corobj<-cor(.z)
cormat<-matrix(corobj,nrow=nrow(corobj))
cormat[upper.tri(cormat)]<-NA
rownames(cormat)<-rownames(corobj)
colnames(cormat)<-colnames(corobj)
sink("Cormat.tex")
print(xtable(cormat))
sink()

###########
#Correlation (econ)
###########
.z<-aic_clean[aic_clean$econ.policies==1,c("pred50_sw","pred90_sw","b.scaled.intgrp","m.scaled.intgrp",
					"rescaled_ideol_code","bioutcome","rep","dem")]
.z<-.z[complete.cases(.z),]
corobj<-cor(.z)
cormat<-matrix(corobj,nrow=nrow(corobj))
cormat[upper.tri(cormat)]<-NA
rownames(cormat)<-rownames(corobj)
colnames(cormat)<-colnames(corobj)
sink("Cormat_Econ.tex")
print(xtable(cormat))
sink()

###########
#Correlation (FP)
###########
.z<-aic_clean[aic_clean$fp.policies==1,c("pred50_sw","pred90_sw","b.scaled.intgrp","m.scaled.intgrp",
					"rescaled_ideol_code","bioutcome","rep","dem")]
.z<-.z[complete.cases(.z),]
corobj<-cor(.z)
cormat<-matrix(corobj,nrow=nrow(corobj))
cormat[upper.tri(cormat)]<-NA
rownames(cormat)<-rownames(corobj)
colnames(cormat)<-colnames(corobj)
sink("Cormat_FP.tex")
print(xtable(cormat))
sink()

###########
#Correlation (NE)
###########
.z<-aic_clean[aic_clean$ne.policies==1,c("pred50_sw","pred90_sw","b.scaled.intgrp","m.scaled.intgrp",
					"rescaled_ideol_code","bioutcome","rep","dem")]
.z<-.z[complete.cases(.z),]
corobj<-cor(.z)
cormat<-matrix(corobj,nrow=nrow(corobj))
cormat[upper.tri(cormat)]<-NA
rownames(cormat)<-rownames(corobj)
colnames(cormat)<-colnames(corobj)
sink("Cormat_NE.tex")
print(xtable(cormat))
sink()

############
#Distributions
############
baseplot<-ggplot(na.omit(vars))+theme_bw()+theme(plot.title=element_text(size=7),
			plot.subtitle=element_text(size=5),
			axis.title=element_text(size=6),
			axis.text=element_text(size=5))

h1<-baseplot+geom_density(aes(pred50_sw))+
		ggtitle("Middle class policy support",
			sub=paste("mean=",as.character(round(mean(aic_clean$pred50_sw),3)),
				";",
				"SD=",as.character(round(sd(aic_clean$pred50_sw),3))
				)
			)
h2<-baseplot+geom_density(aes(pred90_sw))+
		ggtitle("Affluent policy support",
			sub=paste("mean=",as.character(round(mean(aic_clean$pred90_sw),3)),
				";",
				"SD=",as.character(round(sd(aic_clean$pred90_sw),3))
				)
			)
h3<-baseplot+geom_density(aes(b.scaled.intgrp))+
		ggtitle("Business interest group policy support",
			sub=paste("mean=",as.character(round(mean(aic_clean$b.scaled.intgrp),3)),
				";",
				"SD=",as.character(round(sd(aic_clean$b.scaled.intgrp),3))
				)
			)
h4<-baseplot+geom_density(aes(m.scaled.intgrp))+
		ggtitle("Advocacy interest group policy support",
			sub=paste("mean=",as.character(round(mean(aic_clean$m.scaled.intgrp),3)),
				";",
				"SD=",as.character(round(sd(aic_clean$m.scaled.intgrp),3))
				)
			)
h5<-baseplot+geom_histogram(aes(factor(fp.policies)),stat="count")+
		ggtitle("Foreign policies",
			sub=paste("Binary variable. Percent of observations=",as.character(round(mean(aic_clean$fp.policies),3)))
			)
h6<-baseplot+geom_histogram(aes(factor(econ.policies)),stat="count")+
		ggtitle("Economic policies",
			sub=paste("Binary variable. Percent of observations=",as.character(round(mean(aic_clean$econ.policies),3)))
			)
h7<-baseplot+geom_histogram(aes(factor(ne.policies)),stat="count")+
		ggtitle("Social policies",
			sub=paste("Binary variable. Percent of observations=",as.character(round(mean(aic_clean$ne.policies),3)))
			)
h8<-baseplot+geom_histogram(aes(factor(rescaled_ideol_code)),stat='count')+
		ggtitle("Ideological code",
			sub=paste("Ordered categorical variable")
			)+scale_x_discrete(labels=c("Most conservative"," "," ","Neutral"," "," ","Most liberal"))
h9<-baseplot+geom_histogram(aes(factor(bioutcome)),stat="count")+
		ggtitle("Policy passage",
			sub=paste("Binary variable. Percent of observations=",as.character(round(mean(aic_clean$bioutcome),3)))
			)
h10<-baseplot+geom_histogram(aes(factor(rep)),stat="count")+
		ggtitle("Republican policy support",
			sub=paste("Ordered categorical variable")
			)+scale_x_discrete(labels=c("Low support"," ","Neutral"," ","Most support"))
h11<-baseplot+geom_histogram(aes(factor(rep)),stat="count")+
		ggtitle("Democratic policy support",
			sub=paste("Ordered categorical variable")
			)+scale_x_discrete(labels=c("Low support"," ","Neutral"," ","Most support"))

pdf("DescriptiveDistributions.pdf")
grid.arrange(h1,h2,h3,h4,h5,h6,h7,h8,h9,h10,h11,top="Distributions for variables of interest")
dev.off()
